All square roots are periodic when written as continued fractions and can be written in the form:
For example, let us consider √23: √23 = 4 + √23 — 4 = 4 +
If we continue we would get the following expansion: √23 = 4 +
The process can be summarised as follows:
It can be seen that the sequence is repeating. For conciseness, we use the notation √23 = [4;(1,3,1,8)], to indicate that the block (1,3,1,8) repeats indefinitely.
The first ten continued fraction representations of (irrational) square roots are:
√2=[1;(2)], period=1 √3=[1;(1,2)], period=2 √5=[2;(4)], period=1 √6=[2;(2,4)], period=2 √7=[2;(1,1,1,4)], period=4 √8=[2;(1,4)], period=2 √10=[3;(6)], period=1 √11=[3;(3,6)], period=2 √12= [3;(2,6)], period=2 √13=[3;(1,1,1,1,6)], period=5
Exactly four continued fractions, for N ≤ 13, have an odd period.
How many continued fractions for N ≤ 10000 have an odd period?
In [1]:
from math import gcd
def reduce_triple(a, b, c):
d = gcd(a, gcd(b, c))
if c < 0:
d = -d
return (a // d, b // d, c // d)
def cf_period(n):
sq = n ** .5
a, b, c = 0, 1, 1
history = {}
count = 0
while (a,b,c) not in history:
history[a,b,c] = count
count += 1
k = int((a + b * sq) / c)
a, b, c = reduce_triple(c*(a-c*k), -b*c, (a-c*k)**2-b*b*n)
return count - history[a,b,c]
count = total = 0
for n in range(2,10001):
sq = round(n ** .5)
if sq * sq == n:
continue
if cf_period(n) % 2:
count += 1
print(count)
Note: This is old code and it should be revised and documented.